function M = getM(n_seg, n_order, ts)
    M = [];
    for k = 1:n_seg
        M_k = [];
        T = ts(k);
        %#####################################################
        % STEP 1.1: calculate M_k of the k-th segment 
        M_k = [   0,       0,      0,      0,     0,   0, 0, 1;
                  0,       0,      0,      0,     0,   0, 1, 0;
                  0,       0,      0,      0,     0,   2, 0, 0;
                  0,       0,      0,      0,     6,   0, 0, 0;
                T^7,     T^6,    T^5,    T^4,   T^3, T^2, T, 1;
              7*T^6,   6*T^5,  5*T^4,  4*T^3, 3*T^2, 2*T, 1, 0;
             42*T^5,  30*T^4, 20*T^3, 12*T^2,   6*T,   2, 0, 0;
            210*T^4, 120*T^3, 60*T^2,   24*T,     6,   0, 0, 0];
        
        M = blkdiag(M, M_k);
    end
end